Motivated by the development of efficient Monte Carlo methods
for PDE models in molecular dynamics,
we establish a new probabilistic interpretation of a family of divergence form
operators with discontinuous coefficients at the interface
of two open subsets of $\mathbb{R}^d$. This family of operators includes the case of the
linearized Poisson-Boltzmann equation used to
compute the electrostatic free energy of a molecule.
More precisely, we explicitly construct a Markov process whose
infinitesimal generator belongs to this family, as the solution of a SDE
including a non standard local time term related to the interface
of discontinuity. We then prove an extended
Feynman-Kac formula for the Poisson-Boltzmann equation.
This formula allows us to justify
various probabilistic numerical methods to
approximate the free energy of a molecule.
We analyse the convergence rate of these simulation procedures and
numerically compare them on idealized molecules models.